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ABSTRACT 


Two  methods  of  obtaining  sensitivity  data  were  simulated  on  an 
electronic  computer  for  the  purpose  of  comparing  the  accuracy  of  the 
estimates  of  the  parameters  of  an  underlying  cumulative  normal  response 
function.  The  first  method  simulated  the  standard  Bruce ton  procedure 
while  the  second  used  a  modified  binary  search  routine  with  a  portion  of 
the  sample  in  order  to  obtain  maximum  likelihood  estimates  of  the  input 
parameters  for  use  in  a  foilow-on  3ruceton  test. 

The  results  showed  both  methods  ut)  be  effective  in  estimating  the 
mean  but  with  slightly  more  variaDility  in  the  estimates  obtained  by  the 
second  procedure.  Both  methods  underestimated  the  standard  deviation  - 
again  with  more  variability  in  the  estimates  obtained  by  the  second 
procedure.  When  the  prior  parameter  estimates  were  unknown  and  the 
applicable  stimulus  level  bounded,  the  second  method  yielded  estimates 
favorably  comparable  to  those  expected  frcm  the  Bruceton  procedure  with 
suitable  prior  input  estimates . 
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I.  INTRODUCTION 


Frequently  a  statistician  is  faced  with  the  problem  of  determining 
the  level  of  a  stimulus  which  critically  affects  the  performance  of  a 
device.  The  nature  of  the  testing  to  be  discussed  is  such  that  once 
some  positive  level  of  the  stimulus  is  applied  to  the  device  either  a 
response  or  a  non-response  can  be  immediately  observed  ana,  in  either 
case  device  is  altered  so  that  a  bonafide  result  cannot  be  obtained 

fiom  a  second  test.  Tests  of  this  type  are  known  as  sensitivity  tests. 

One  of  the  many  problems  besetting  those  involved  in  explosives 
research  is  that  of  providing  measures  and  specifying  rules  to  provide 
for  the  safe  handling  and  transportation  of  explosives.  Many  different 
types  of  sensitivity  testing  apparatus  have  been  developed  for  laboratory 
use,  the  most  common  being  those  that  subject  some  quantity  of  explosive 
to  the  impact  load  of  a  falling  drop-weight  from  some  controllable  height. 
At  least  as  late  as  October  1965  there  remained  two  important  physical 
problems  to  be  solved;  namely,  that  of  establishing  a  measure  of  stimulus 
not  highly  apparatus -dependent  and  then  that  of  translation  of  these 
results  to  safe  handling  rules  [1].  These  problems  are  not  addressed  in 
this  paper  but  should  be  kept  in  mind  when  considering  the  overall  problem. 

In  the  early  1940's,  a  technique  for  obtaining  sensitivity  data  was 
developed  and  used  in  explosives  research  at  the  Explosives  Research 
Laboratory,  Bruceton,  Pennsylvania  which  has  come  to  be  called  synony¬ 
mously,  the  Bruceton,  Staircase,  or  "Up  and  Down’’  Method. 

The  aim  of  this  method  of  testing  is  to  increase  the  accuracy  with 
which  certain  critical  values  of  the  stimulus  may  be  estimated,  notably 
the  median  (or  mean)  and  standard  deviation.  The  accuracy  of  the  method 

Preceding  page  blank 
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depends  in  part  on  the  stimulus  level  at  which  the  first  item  is  tested 
and  the  interval  spacing  for  subsequent  levels  of  testing  [2]. 

When  the  stimulus  levels  mentioned  above  cannot  be  determined  prior 
to  testing  or  when  little  confidence  is  placed  on  the  available  estimates, 
a  preliminary  (or  search)  phase  of  testing  may  be  desirable  to  obtain 
maximum  likelihood  estimates  prior  to  employing  the  Bruceton  Method  with 
the  remainder  of  the  sample.  A  procedure  to  do  this  is  offered  as  an 
alternative  method. 

The  comparative  accuracies  of  the  two  techniques  were  examined 
through  the  use  of  simulation  conducted  on  a  high-speed  electronic 
computer.  All  parameters  and  estimates  considered  as  inputs  to  the 
simulation  were  kept  within  ranges  for  which  the  Bruceton  Method  is 
considered  to  yield  accurate  results  [2], 


II.  THS  MODEL 

Let  x  be  an  applied  stimulus  level  (xe[o,®))  and  y  =  y(x)  be  the 
associated  response  (ys<^o,l^  where  "o"  denotes  no  response  and  “1" 
denotes  response).  At  any  given  stimulus  level  consider  y  to  be  the 
realization  of  a  Bernoulli  random  variable,  Y,  with  response 
probability 

p(x)  =  Prob  (Y  -  ljx) 

The  function  p(x)  is  called  the  response  function  and  is  further 
specified  as 

p(x)  »  0  xe[o,a] 

0  <  p(x)  <  1  xe(a,b) 

and  p(x)  «  1  xe[b,<=) 

The  intervals  [o,a],  (a.b)  and  [b,®)  are  called  the  zero-response 
region,  the  mixed-response  region,  and  the  one-response  region 
respectively.  It  is  assumed  that  p(x)  is  a  saonotonely  increasing 
function  for  stimulus  values  in  the  mixed-response  region.  Thus, 
p(x)  can  be  considered  as  the  cumulative  distribution  function  for  a 
random  variable  X  such  that 

p(x)  =>  Prob  <X  <  x).  [3] 

In  this  context  the  random  variable  X  can.  be  interpreted  as  a  thres¬ 
hold  stimulus  level,  thus 

Prob  (Y  =  ljx)  =  Prob  (X  <  x)  =  p(x) 

and  Prob  (Y  =  ojx)  **  Prob  (X  >  x)  «*  1  -  p(x).  [33 

2 

It  is  assumed  the  X  is  distributed  Normal  (,j,0  );  that  is 

p(x)  »  {p(xJn,o  ) 

o 

where  cpCxjjijO*’)  represents  the  cumulative  normal  distribution  with  mean 
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and  variance  e  .  In  particular 


III.  TESTING  METHODS 


A.  BRUCETON  METHOD 
1.  Description 

Based  on  intuition  or  past  experiments,  the  experimenter  selects 
a  priori  estimates  of  |i  and  a.  Call  these  estimates  and  Oj.  and  let 


The  experimenter  tests  the  first  item  at  or  near  (ij.  If  there 
is  a  response  the  second  item  is  tested  at  a  level  d  units  below  |i^, 
otherwise  the.  second  item  is  tested  at  a  level  d  units  above  ;i  .  In 
the  same  manner,  each  of  the  remaining  items  is  tested  at  a  level  d 
units  above  or  below  the  previous  test  level  according  as  there  was  not 
or  there  was  a  response  observed  for  the  previous  test.  Thus  the 
sample  is  concentrated  about  the  mean  and  one  would  expect  nearly  equal 
numbers  of  responses  and  non-responses.  In  fact,  the  number  of  non¬ 
responses  at  any  level  will  not  differ  by  more  than  one  from  the 
number  of  responses  at  the  next  higher  level  [2]. 

Let  N  denote  the  total  number  of  observations  of  the  less 
frequent  event  and  nQ,n^,n2 , • • *n^  denote  the  frequencies  of  this  event 
at  each  level  where  n^  corresponds  to  the  lowest  level  and  n^  the 
highest  level  at  which  the  less  frequent  event  occurs. 

The  final  estimates  of  }i  and  o  are  based  on  the  first  two 
moments  of  the  stimulus  levels.  Since  the  intervals  are  equally 
spaced,  these  moments  can  be  computed  in  terms  of  the  sums 

A  =»  S  i  n. 

i 


and 


B  o  S  i 
i 


n. 
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Let  ji  be  the  estimate  of  p,  by  this  method.  Then 


£  -  *’ +  <*  d  ±  D 

where  x*  represents  the  lowest  level  at  which  the  less  frequent  event 
occurs  [2].  The  plus  sign  is  used  when  the  analysis  is  based  on  non¬ 
responses,  and  the  minus  sign  when  it  is  based  on  responses  [2]. 

2  2 

If  (NB-A  )/N  >  .3  the  sample  standard  deviation  is 

2 

s  =»  1.620  d  +  .029) 

Otherwise,  a  more  elaborate  calculation  must  be  employed  and  is 
described  in  Ref.  2. 

To  obtain  confidence  intervals,  estimates  of  the  standard 
deviations  of  the  sample  mean  and  sample  standard  deviation,  say  s^ 
and  s  ^  respectively,  are  given  by 

b 


s 


ta 


Gs 


and 

Hs 

e  «  — 

8  ,/ii 

where  the  factors  G  and  H  are  dependent  on  the  ratio  ^  and  the 
position  of  the  mean  relative  to  the  testing  levels.  Plots  of  these 
factors  are  available  in  Ref.  1. 

2.  Discussion 

Only  rarely  is  the  threshold  stimulus  Z  normally  distributed. 
It  is  usually  the  case  that  some  scale  transformation  of  Z,  say  X,  is 
made  so  that  X  is  normally  distributed  in  the  vicinity  of  the  :nean. 
This  transformation  is  done  prior  to  testing  to  determine  and 
Only  after  all  analysis  is  completed  are  the  values  scaled  back  to  the 
original  stimulus  measure  [2}. 
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The  size  of  the  sample  is  critical  to  the  accuracy  of  the 
estimation.  Note  that  at  most  only  half  of  the  sample  is  used  in  the 
analysis  so  that,  for  example,  if  thirty  items  are  tested  the  maximum 
possible  value  of  N  is  fifteen.  The  analysis  is  based  on  large  sample 
theory  which  in  the  case  mentioned  would  be  applied  to  a  sample  of  size 
fifteen  [2]  [4]. 

Unless  normality  of  the  variate  is  assured  this  method  does 
not  yield  accurate  results  for  the  small  and  large  percentage  points. 
This  is  unfortunate  since  in  most  applications  one  would  be  more 
interested  in  a  small  percentage  point  as  a  measure  of  safety  and  a 
large  percentage  point  as  a  measure  of  reliability.  At  any  rate,  an 
estimate  of  a  percentage  point  j  is 

j  »  Jji  +  ks 

where  k  is  chosen  from  tables  of  the  standard  normal  deviate  to  give 
the  desired  percentage  t2].  One  could  then  conduct  tests  in  the 
vicinity  of  this  value  to  refine  the  estimate. 

B.  BRHCETON  METHOD  PRECEDED  BY  SEARCH 
1.  Description 

In  the  event  that  a  priori  estimates  of  pi  and  o  are  not 
available  some  economic  method  of  attaining  these  estimates  is  desired. 
A  method  proposed  and  described  below  is  a  modified  binary  search 
technique . 

Again,  the  assumption  is  that  the  threshold  stimulus  (or 
some  transformation  of  it)  is  normally  distributed  and  p(x)  can  be 
represented  by  a  cumulative  normal  distribution. 

As  noted  from  the  model 

Prob  (Y  13  ojx  <  a)  1 
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and 


Prob  (Y  =  1 1  x  >  b)  «  1. 

The  first  step  in  the  procedure,  then,  is  to  select  values  for  a  and  b. 
(In  the  case  of  complete  uncertainty  these  could  be  the  limiting  values 
of  the  testing  apparatus)  and  commence  the  binary  search  starting  at 

x  =»  (a  +  b)/2. 

If  p(x)  were  a  step  function,  repetition  of  this  method  would  locate 
the  step  in  an  interval  of  any  desired  length-  In  general,  however, 
the  mixed-response  region  has  non-zero  width  and  a  non-response  would 
merely  indicate  that  the  applied  stimulus  is  in  the  mixed  response 

region  or  below  while  a  response  would  indicate  that  it  was  in  t.he 

mixed  response  region  or  above. 

If  a  test  at  yields  a  response  and  a  test  at  x^  yields  a 

non-response  while  x^  <  x^  it  is  certain  that  both  x^  and  x^  are  in  the 

mixed  response  region.  This  condition  is  called  a  response  inversion 
and  is  the  basic  indicator  for  the  modified  binary  search  technique. 

The  description  of  the  procedure  is  best  followed  by  referring  to 
Figures  1  through  4. 

Sequence  S*  is  a  cyclic  one  indicating  that  a  reduction  in  step 
size  should  be  taken.  Test  levels  are  selected  attempting  to  reproduce 
this  sequence.  Failure  to  do  this  results  in  the  basic  inversion 
sequence  Sq.  Tests  are  then  made  at  the  end  of  this  sequence  to  result 
in  one  of  three  terminal  situations  S^,  S^,,  or  S^.  In  the  event  the 
mixed  response  region  is  relatively  narrow  and  near  a  or  b,  several 
binary  reductions  may  be  necessary  to  reproduce  S*  or  one  of  the 
terminal  situations.  These  circumstances  are  represented  by  and 
Sjj  [3]. 
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STARTING  SEQUENCE 
Figure  1 
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Maximum  likelihood  estimates  of  ^  and  a  are  available  for 
sequences  S^,  S^,  and  and  developed  as  described  b 2 low  [3] . 

2.  Discussion 

It  is  assumed  that  all  trials  are  independent.  Thus  the 
probability  of  the  sequence  is 

Prob  (S1)  -  Prob  (Y j=0 ,  Y2=0,  Y3=l,  Y^O,  Y^l,  Y^lIx^.X^X^ 


) 


6 

=  tt  Prob  (Y.  =*  y.  |x.  > 
i=l  x  1  1 


where 


and 


Prob  (Y^  53  y1|xi)  53  cpO^)  if  y.  “1 

=  1  -  cp(xi)  if  y£  =  0 


cpC^) 


»  Prob  (X.  <  x.) 

1  —  X 


1  (X-U)2 

2  2 


dx. 


Maximum  likelihood  estimates  for  u  and  a  can  then  be  established 
using  standard  normal  tables  for  each  of  the  terminal  situations. 

These  estimates  are  indicated  on  Figure  5. 
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RESPONSE  SEQUENCES  AND  PARAMETER  ESTIMATES 
Figure  5 
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IV.  SIMULATION 


A.  DESCRIPTION 

All  simulated  experiments  were  conducted  on  an  IBM  360/67  computer 
using  the  FORTRAN  IV  programming  language.  The  basic  program  is 
attached.  The  response  function  p(x)  used  was  cumulative  normal  with 
[1  °  30  and  o  =  3. 

The  sample  size  was  kept  at  seventy  for  each  experiment  to  provide 
some  assurance  that  the  analytical  sample  would  be  suitable  for  large 
sample  ana lys is . 

The  basic  test  procedure  was  to  draw  a  random  number  on  the  unit 
interval  and  compare  this  to  F(x),  a  function  of  a  standard  normal 
variate  specified  as 


and 


where 


if  x  <  0 


if  x  >  0 


erf (v) 


(The  function  subprogram  erf  is  an  IBM-supplied  subprogram.)  If  the 

random  number  was  less  than  or  equal  to  F(x^)  then  a  response  was 
th 

counted  for  the  i  level;  otherwise  a  non-response  was  counted. 

Six  different  cases  were  tested  using  the  straight  Bruceton  pro¬ 
cedure  (METHOD  1)  with  two  different  input  estimates  of  |i  and  three 
different  input  estimates  of  a.  Case  1  considered  exact  estimates; 
i.e.,  VJLj  «  M  and  a ^  a.  Case  2  considered  =  [_i-6  and  Oj  **  o* 
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Cases  3  and  4  considered  a  p  and  a^  -  a/2,  2a  respectively  while 
Cases  5  and  6  repeated  Cases  3  and  4  except  p^.  =  p-6.  For  each  of  the 
six  cases  1000  experiments  were  conducted  each  utilizing  a  different 
sequence  of  random  numbers . 

The  search  procedure  (METHOD  2)  was  then  incorporated  into  each 
of  the  above  six  cases  using  the  a  prior  estimates,  p^  and  a_,  to 
determine  estimates  for  stimulus  levels  a  and  b  and  thereby  the  sizo 
of  the  binary  reduction  as  indicated  in  Figure  1.  The  program  then 
followed  the  flow  shown  in  Figures  1  through  4  until  either  a  terminal 
sequence  was  reached  or  the  search  was  arbitrarily  terminated  as  dis¬ 
cussed  in  subparagraph  C  below.  The  Bruceton  procedure  was  then  used 
until  the  sample  was  exhausted. 

The  final  case.  Case  7,  indicated  complete  lack  of  knowledge  of  p 
and  0  but  considered  the  upper  and  lower  stimulus  level  limits  of  the 
test  apparatus  to  be  100  and  0  respectively. 

B.  MEASURES  OF  EFFECTIVENESS 

At  the  completion  of  all  experiments  for  each  case,  several 
measures  were  obtained  for  comparison.  First,  average  values  of  the 
parameters  were  determined  to  be 

p  »  E  p./N 
i 

and 

a  =  S  a./N 
i 

where  p^  and  are  the  a  posteriori  estimates  of  p  and  a  for  the  i*^ 
experiment  and  n,  the  number  of  experiments  used.  Next,  as  measures  of 
variability 
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s2a  "=  E  (|J.  -  n)2/n-l 

i  i 


and 


2 

s  a  »  E 
0  i 


o)/n-l 


were  calculated.  In  addition,  the  program  listed  the  maximum  and 
minimum  estimates  of  both  j_i  and  g. 


C.  DISCUSSION 

In  Chapter  III  it  was  noted  that  sequences  S*,  StT,  and  ST  are 

U  Li 

cyclic.  In  order  to  simplify  the  program  it  was  necessary  to 
artificially  terminate  these  situations  at  some  point  and  calculate 
the  input  values  for  the  Bruceton  test.  The  estimate  of  (a  used  was 

^  *  <X1  +  *2>12 

where  x^  and  x^  are  adjacent  testing  levels  and  x^  >  x^  with  y^  =  0 
and  y2  -  I.  The  estimate  of  a  used  was 

as  "  (x2  "  xl)/2 
for  Cases  1  through  6  and 

as  =  (x2  '  xl)/6 

for  Case  7.  The  former  estimate  of  o  was  chosen  arbitrarily  while  the 
latter  estimate  was  based  on  the  estimate  of  the  mixed  response  region 
being  6a.  While  the  number  of  terminations  of  this  type  was  insignifi¬ 
cant  for  the  first  six  search  cases,  in  the  final  case  over  600  experi¬ 
ments  were  thus  terminated  requiring  the  program  to  be  expanded  to 
permit  more  recycling.  The  point  is  that  the  artificial  termination 
does  not  represent  the  search  procedure.  This  problem  would  not  arise 
in  field  experimentation  until  either  the  sample  was  exhausted  or  the 
step  size  reduction  of  stimulus  level  indicated  was  too  narrow  to  be 
measur  j  or  controlled  by  the  test  apparatus. 
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Also  in  the  interest  of  program  simplification  those  experiments 


for  which 


<  .3 


were  not  used  for  analyses.  This  limitation  invalidated  the  measures 
of  effectiveness  for  the  Bruceton  cases  where  o^.  -  2cr. 


D.  RESULTS 

The  results  of  the  simulation  are  listed  in  Table  I.  It  is 
questionable  that  the  measures  listed  under  Method  1  are  valid  for 
Cases  4  and  6  in  that  only  .381  and  .393  of  the  possible  experiments 
were  used.  These  two  cases  and  Case  4  under  Method  2  (where  .661  of 
the  possible  experiments  were  used)  are  the  only  ones  for  which 
a  >  a. 

In  general  the  extreme  estimates  are  more  widely  separated  and  the 
variability  of  a  is  greater  in  Method  2. 

Estimates  of  (i  range  from  27.8823  to  31.7647  for  Method  1  and 
27.937  to  31.91  for  Method  2. 

Estimates  of  o  range  from  .8741  to  6.5027  for  Method  1  and  .3498 
to  9.8328  for  Method  2. 

The  lowest  average  (i,  29.9113,  was  obtained  under  Method  1,  Case  5, 
while  the  highest  average  0,  30.1175,  was  obtained  under  Method  2, 

Case  3. 

The  lowest  average  a,  2.3748,  was  obtained  under  Method  2,  Case  5, 
while  the  highest  average  O,  2.9474,  was  obtained  under  Method  1, 

Case  5.  (Case  6  is  not  counted  under  Method  1  nor  is  Case  4  under 
both  methods.) 
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TABLE  OF  EXPERIMENTAL  RESULTS 


METHOD  1 

METHOD  2 

A 

A 

A 

1  A 

0 

u 

a 

29.9641 


31.6765 


28.3235 


2 .9040 


5.8249 


31.4571 


28.0286 


.2574 


3.5424 


6.0266 


.4639 


CASE 

5 

^24 

AVE 

Oj-  3  1.5 

MAX 

b  =  30  1 

I  VAR 

CASE 

6 

1.4773 


.2220 


2  9474 


5.9257 


29.9493 


31.4118 


27.8823 


.2639 


3.5438 


=  0 


=  100 


30.0317 


31.6875 


28.1976 


29.9363 


31.4063 


20.1961 


30.0247 


28.2552 


30.0123 


31.8229 


27.9541 


2.8609 


5.9343 


1.1241 


.5831 


2.8819 


9.1369 


.9512 


.6336 


Hj-30 

AVE 

30.0295 

2.7216 

30.1175 

2.8615 

a-j.  *  1.5 

MAX 

31.6071 

6.1197 

31.9100 

7.7997 

a  »  24 

MIN 

28.4118 

.8741 

28.5950 

.8697 

b  «  36 

VAR 

.2046 

.7409 

.2693 

.9081 

29.9750 

3.0721 

31.6875 

6.3569 

27.9370 

1.6170 

.2619 

.4522 

2.3748 


7.9507 


.3498 


1.4889 


2.3398 


6.4201 


9.8328 


.5082 
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V.  CONCLUSIONS  AND  RECOMMENDATIONS 


A  CONCLUSIONS 

1.  Estimation  of  the  Mean 

Both  methods  estimate  the  mean  effectively. 

2.  Estimation  of  the  Standard  Deviation 

Both  methods  tend  to  under-estimate  the  standard  deviation  with 
no  predictable  bias  and  are  therefore  unsuitable  for  use  in  safety  or 
reliability  statements.  This  conclusion  agrees  with  the  findings  of 
Hampton  [4]  as  it  pertains  to  the  Bruceton  Method. 

3.  Extension  of  the  Search  Phase  for  the  Starting  Sequence 

Termination  of  the  search  phase  with  sequence  S^  in  the  starting 

sequence  (see  Figure  1)  may  yield  estimates  of  a  greater  than  twice 
the  actual  value.  To  avoid  this  it  is  advisable  to  extend  the  search 
phase  as  described  in  Ref.  3. 

4.  Use  of  Search  Technique 

The  search  procedure  should  be  used  in  thosa  cases  where  there 
is  not  independent  evidence  that  the  estimate  of  a  is  within  the  range 
for  which  tha  Bruceton  Method  is  recommended  (i.e.,  a/2  <  <  2a). 

B.  RECOMMENDATIONS 

Further  testing  of  Method  2  is  recommended  under  the  circumstances 
listed  below. 

1.  Reduction  of  Sample  Size 

It  would  be  of  interest  to  reduce  the  sample  size  to  the  point 
vhere  the  effective  sample  is  small,  say  15,  and  compare  the  Bruceton 
procedure  with  the  search  procedure  using  the  entire  sample  for  the 
search . 
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2.  Random  Selection  of  Response  Function  Parameters 

A  more  valid  test  of  both  methods  would  be  achieved  by 
randomally  selecting  values  of  ji  and  a  over  some  range  and  using 
the  on-line  computer  facility  to  conduct  the  simulation. 
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COMPUTER  PROGRAM 


THIS  PROGRAM  SIMULATES  SENSITIVITY  TESTING  BY  BOTH  THE  BRIfC- 
ETON  METHOD  (WHEN  IANY=0)  AND  THE  RRUCETON  METHOD  PRECEDED 
BY  THE  MODIFIED  BINARY  SEARCH  (WHEN  I ANY= I )  .  THE  UNDERLYING 
RESPONSE  PUNT  ION  IS  CUMULATIVE  NORMAL  <30,3).THF  im^UT  EST¬ 
IMATES  OF  THE  MEAN  AND  THE  STANDARD  DEVIATION  ARF  CALLED 
EXMU  AND  EXSIG  respectively. 

THE  PRINCIPLE  VAP TABLE  NAMES  ARE  AS  FOLLOWS... 

AACT  IS  THE  STIMULUS  VALUE  AT  THE  UPPER  LIMIT  OF  THE 
MIXED  RESPONSE  REGION. 

8ACT  IS  THE  STIMULUS  VALUE  AT  THE  LOWER  LIMIT  OF  THE 
MIXEL  RESPONSE  REGION. 

A  AND  B  ARE  ESTIMATES  OF  AACT  AND  RACT  RESPECTIVELY. 

XU)  IS  THE  STIMULUS  LEVEL  OF  THE  JTH ,  STIMULUS. 

IXOIJ)  IS  THE  CUMULATIVE  COUNT  OF  NON-RESPONS FS  AT  X(J). 
IXX(J)  IS  THE  CUMULATIVE  COUNT  OF  RESPONSES  AT  XU) 

IS  IS  THE  SAMPLE  SI7E. 

NU  IS  THE  ENTRY  NUMBER  FOR  THE  RANDOM  NUMBER  GENERATOR, 
UNIF. 

N  COUNTS  THE  NUMBER  OF  EXPERIMENTS. 

RN  IS  THE  RANDOM  NUMBER  ON  (0,1)  RETURNED  BY  UNIF. 

FOFX  IS  THE  VALUE  OF  THE  RESPONSE  FUNCTION  RETURNED  BY 
SUBPROGRAMS  XNCOF  AND  SNCDF. 

I  SUMO  IS  THE  TOTAL  NUMBER  OF  NON-RESPONSES  FOR  ONE  EXPER¬ 
IMENT. 

ISUMX  IS  THE  TOTAL  NUMBER  OF  RESPONSES  FOR  ONE  EXPERIMENT 
NT  IS  THE  MINIMUM  Oc  TSUMO  AND  ISUMX. 

NSJJ )  IS  THE  FREQUENCY  OF  THE  LESS  FREQUENT  EVENT  AT  X(J) 
NG(J)  REARRANGES  NS(J)  SO  THAT  NGfl)=NSm  WHERE  X(I)  IS 
THE  LOWEST  STIMULUS  LEVEL  AT  WHICH  THE  LESS  FREQUENT 
EVENT  OCCURS. 

APU)  IS  USED  Tn  CALCULATE  THE  FIRST  MOMENT  ,SUMAP. 

8RU)  IS  USED  TO  CALCULATE  THE  SECOND  MOMENT »  SUMBR . 

YPRIME  IS  THE  LOWEST  LEVEL  AT  WHICH  THE  LESS  FREQUENT 
EVENT  OCCURS 

XMUEST  IS  THE  FINAL  ESTIMATE  0=  THF  TRUE  MFAN •  XMtJ. 

DEVEST  IS  THE  FINAL  ESTIMATE  OF  THE  TRUE  STANDARD  DEVI AT- 
ION,XSIG. 

EMUCJ)  IS  THF  DIFFERENCE  OF  XM'JEST  AND  XMtJ. 

EOEVU)  IS  THE  DIFFERENCE  OF  DEVEST  AND  XSIG. 

SAMAVM  AND  SAMAVD  ARE  THE  SAMPLE  AVERAGF  ERRORS  OF  XMUEST 
ANO  DEVEST  RESPECTI VFLY. 

SAMSQM  AND  SAMSQD  ARE  THE  AVERAGE  mcaN  SQUARE  eRRORS  OF 
XMUEST  AND  D6VFST  RESPECTIVELY. 

NOGO  IS  THE  NUMBER  OF  EXPERIMENTS  NOT  USED  IN  THE  FINAL 
ANALYSIS 

DIMENSION  ARRAYS  AND  FORMAT 
SIMULATE  RRUCETON  FIRST  THEN  SEARCH 

IANY-0 

69  TH ING=0. 

INITIALIZE  INTERNAL  AND  OUTPUT  VARIABLES 

SET  EXPERIMENT  COUNTER  ,  SAMPLE  SIZE  COUNTER,  AND  NU- 

LCOUNT =t 
NU=1237l 
103  N=1 

IF( IANY.EQ.O)  N8R=0 
I F ( IANY.EQ.l)  NBR=1 

SET  INPUT  VARIABLES 

XMU=30, 

XSIG=3. 

A-0 . 

B=1 00. 

EXMU=50. 
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EXSIG=12.5 
A=  EXMU-IQ*EX$TG 
B=  FXMU*IQ*EXSIG 
Xl  = { A+B ) /2 • 

10=4 
1NC  =  0 


PROVIDE  BRANCH  TO  STANDARD  BRUCETON 
IF(NBR.EO.O)  GO  TO  33 

CONDUCT  SEARCH 


CALL  UNIfMRN,NU) 
FOFX=XNCDF { XI  *  XHU» 
IFlRN«.LE*FOPX)  GO 
X2»(  B  +  XD/2. 
NBR=NBR+I 
CALL  UN  IP ( RN  «MU ) 
FOFX  =  XNCOF(  X?,XM(J, 
IFt&NoLE. F0C  X )  GO 
X3=(R+X2)/2. 
NBR=NBR H 
CALL  UN  I  F( RN , NU ) 
FOFX=XNCDF( X3f XMUt 
IF(RN.LE.FQPX)  GO 
X4  =  ( B+X3  j /2i 
NBR=NBR  +  'l 
CALL  UNIF(RN,NU I 
FOFX=XNCDF { X4 , XMU* 
IFtRN.LE. FOFXI  GO 
X5=I R+X4J/2. 
NBR=NBR+I 
CALL  UNTF(RN*NII) 
FOFX=XNCOF{  X5,XMIJ, 
IF(RN.LE.FOFX)  G'. 
EXMU={ B+X5J /2. 
£XSIG=  (  X.5-X4)  /?. , 
EXSI G=2**FXSIG 
EXSIG=EXSIG/6c 
GQ  TO  7000 
1313  EXMU  =  tX5+XM/2. 
6XSIG=( X5-X4I/2. 
EXSIG=2.*EXS!G 
EXSIG=EXSIG/6. 

GO  TO  7C00 
9063  X5=( X3+X2) /2. 


XSIG) 

TO  9300 


XSIG) 

TO  9250 


XSIG) 

TO  9125 


XSIG) 

TO  9063 


XSIG) 

TO  1313 


NBR=NBR+1 
CALL  UNIFI RN? NU) 
FDFX=XNCOFf  X5,XMU,X$IG) 
IF(RN.LE.FQFX)  GO  TO  1314 
X6={  X3FX4J/2. 

NBR=NBR+1 
CALL  UNIF{RN,MU) 
FOFX=XNCDF  (  X6  *  XM1J *  XSI G  ) 
IF(RN.LE.FOFX)  GO  TO  1316 
EXMU=<  X6+X4)/2. 
EXSIG=(X6-X3)/2. 
EXSIG=2.*EXSIG 
EX$IG=EXSIG/6* 

GO  TO  7000 
1316  FXMU=( X6+X3) t 2. 
E.XSIG=(X6-X3)/?, 

EXSI G=2 .+EXSIG 
EXSIG=EXS IG/6 . 

GO  TO  7000 
1314  X6=2.*X2-X5 
NBR=NBR+1 
CALI  UN Ic ( RN  *  NU ) 
F0FX=XNCDF(X6vXwU»XSIG) 
IFIPN.LE. EOF X *  GO  TO  1315 


XB=X5 
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»  XSIG) 

TO  909*- 


tXSIG) 

TO  9047 


»XSTG ) 

TO  9024 


♦  XSir, ) 

TO  9012 


DELX=X5-X2 
eXMU=XB+DELX/2. 

EXSI G= 1 .  3*DELX 
GO  TO  7000 
1315  XB  =  Xfe 

DELX=X2-X6 
EXMU=XB-0ELX/4» 
EXSIG=6*0ELX 
GO  TO  7000 
9125  X4=(X2+Xl)/2. 

NBR=NBR+1 
CALL  UNIF(RN,NU> 
FOFX=XNCDF(  X*-rX«U 
inPN.LE.FOFX)  GO 
X5=(  X2+X3J/2. 

NBR=NRR+1 
CALL  UN  IF  (RN,  Nil) 
FOFX=XNCOF( X5,XMU 
IF(RN.LE.FOFX)  GO 
X5=(BfX3) /2. 

NBR«NBR+1 
CALL  UNIF(RN,NU) 
FOFX=XNCOF( X6,XMU 
IFIRN.LF.FOFX)  GO 
X7=?.*B-X6 
NBR=NBIJ  <1 
CALL  UNJF(RN,NU) 
F0FX-XNCDF(X7,XMIJ 
IF(PN.LF.  F0P  X )  GO 
XB=X7 
0ELX=X8-R 
EXMU=XB+0El.X/4« 

FxsrG=6.*oeix 

GO  TO  7000 
9012  XB=X3 

DEL  X  =  X6— XR 
EXMIj=XR+0ELX/2. 

EXSI G= /  3+QELX 
„  GO  TO  7000 
9024  FXMU-I X3+X5J/2. 
EXSIG=(X3-X5)/2. 
EXSIG=2.*FXSIG 
EXSIG=FXS!G/6. 

GO  TO  7000 
9047  X6=(  X4+X2) /2. 

NBR=NBR+1 
CALL  UNIFIRN, 
F0FX=XNC0F(X6 
IF(RN.LE.FOFX)  GO 
EXMU=( X5+X2) /2. 

EXSIG= { X5-X2 ) /2 . 

EXSI G=  2  XS IG 
EXSIG=EXSIG/6. 

GO  TO  7000 
9C1 1  X?=2.*X4-X6 
NBR=NBR*1 
CALL  UN  I F { RN  tNU  J 

fgf&W^l^.XMU.XSIGI 

IFCRN.LE.FOFX)  GO  TO  9010 
Xo  =  X6 

0ELX=X2-XB 
FXMU=XR*OEl.X  f2» 
rXSIG=1.3*neLX 
„  GO  TO  7000 
9010  XB=X7 

DELX=X4~X7 
EXMU=XB“0ELX/4, 
EXSIG=6*0Elx 
GO  TO  7000 
9094  X5=2*XI-X4 
NBR=N8R+ 1 
CALL  UNIF{RN,NU> 


HJ) 

>  t  XMU 


tXSIG) 

TO  9011 
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FOF  X  =  XNCOF { X5  t  XMU 
IFIRN.LE.FOFX)  GG 
XB  =  X4 

DEL X=X2-X4 
EXMIJ=XR+DELX/2. 
EXSIG=l.3*OELX 
GO  TO  7000 


9003  X8=X5 

DELX=Xl-X5 
EXMU=XR-DFLX/4. 
EX$IG=6.*DELX 
GO  TO  7000 
9250  X3= ( A+Xl ) /2. 
NBR=NBR+1 
CALL  UN  IF  ( RN  *NU ) 
FOFX  =  XNCOF(X3»XMi| 
IFIRN.LE.FOPX)  GO 
X4  =  (  XI  +X2)/2. 
NBR=NBR+1 
CALL  UN!F(RN«NU) 
FOF  X=XMC0F(X4»XMUi 
IF(RN.LE.FOFX)  GO 
X5=(B+X2)/2. 
NBR=NBR+l 
CALL  UNIF { RNtNU ) 
c0FX=XNCDF{X5,XMU, 
IFIRN.LE. FOFX)  GO 
X6=2*B-X5 
NBR-NBR+1 


CALL  UN  IF  { RN  »N!J ) 
F0FX'=XNCDF(X6,XMU. 
IF(RN.Lc.FOFX)  GO 
XB- X6 
DELX=X6-R 
EXMU=XB+D6LX/4. 
EXSIG=6*DELX 
GO  TO  7000 
9006  XR=X2 

DELX=X5-XR 
EXMU=XB+DELX/2. 
EXSI G= 1  • 3*DEL  X 
GO  TO  7000 
9005  EXMU=( X2+X4) /?. 
EXSIG=(X2~X4)/2. 
EX$iG=2.*EXSIG 
EXSIG=EXS I G/6. 

GO  TG  7000 
9004  X5=( Xi+X3)/2. 
NBR=NBR+l 
CALL  UNIF( PN.NU) 
F0FX  =  XNCDF(X5,X'Mlf 
IF(RN.LF,FQPXJ  GO 
EXMU  =  {  XHXAI  /2, 
EXSIG= I XA-X1 )  /  2 . 
EXSIG=2.*EXSIG 
EXSI G=F  XS IG/6. 

GO  TO  7000 
5555  X6  =  2.*X3-X5 
NBR=NBR+1 
CALL  UNIF(R*UN<J) 
FOF  X=XNCDF ( X6  »  XMU t 
IFIPN.l EeFOFX)  GO 
X8  =  X5 


DELX=Xl-X^ 
EXMU=XB+DELX/2. 
FXSIG=1.3*PELX 
GO  TO  7000 
9007  XR=X6 

DELX=X3-XB 
EXMU=XB-DFLX/4. 
EXSTG=6*DELX 
GO  TO  7000 


»  X  $  I G  ) 

TO  9003 


XSIG) 

TO  93  75 


XSIG  I 
TO  9004 


XSIG) 

TO  9005 


XSTG) 

TO  9006 


XSIG) 

TO  5555 


XSIG) 

TO  9007 
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9375  X4=?.*A-X3 
NBR=N8R+1 

CALL  UNIFCRN/NU) 
FOFX=XNCDF( X4,XMU 
IF(RN.LE.FOFX)  G  n 
XB=X3 

DEL  X  =  X1~XR 
EXMU=XR+OELX/2. 
EXS  f G= 1 . 3*DE  LX 
GO  TO  7000 

9376  X8=X4 

DEL  X=A-XB 
£XMU=XR-D£LX/4. 
EXSI0=6*D6LX 
GO  TO  7000 
9500  X2  =  <  A  +  XD/2. 
NBR=N8R+1 
CALL  UNIF(  RN  »Ntl } 
F0FX=XNC0F(X2,XMU 
IF(RN.LE.FOFX)  GO 
X3={ X1+8J/2. 

NBR  =  Nr>R  +1 

CALL  UN  IF  (RN, Nil) 

FOFX=XNCCF ( X3 ♦ XMU 

IFC RN.LE.FOFX)  GO 

X4=2.*R-X3 

NBR=NBR+1 

CALL  UMIF(RN,NU) 

FOFX=XNCDF( X4  tXMU 

IFCRN.LE.FQFX)  GO 

XB=X4 

DELX=X4~XR 
EXMU=XB-s-0ELX/4. 
EXSIG=6*0ELX 
m  GO  TO  7000 
5554  XB=X1 

DELX=X3-XB 
EXMU=XB+0ELX/2. 
EXSiG=1.3*0ELX 
GO  TO  7000 

5556  X4=(Xl*X2)/2. 
NBR-NBR+1 

CALL  UNIF{ RN,NU) 
F0FX=XNCDF(X4,XMl? 
JFCRN.LE.FCFX)  GO 
X5=( X1+X3) /2. 
NBR=NBR+1 
CALL  IJN!FCRN9NU) 
FOFX=XNCDF  {  X5 » XMIJ 
IFCRN.LE.FO^x)  go 
Xf-?.*X3-X5 

nbr=;:br+i 

CALL  UNIFCRM,NU) 
FQFX=XNCOr  C  X6»  XMU 
IFCRN.LE.FOFX)  go 
XP.=X  6 

DELX=XB-X3 
EXMU=XR4-0FLX/4. 
EXSIG=6*DELX 
GO  TO  7000 
:>559  X8=  XI 

D5LX=X5-X1 
EXMU=XB+DELX/?. 
EXSI G= I •  3*DELX 
GO  TO  7000 
5553  EXMU=C  X4^X1 ) /2. 
EXS!G=CXl-X4)/2. 
EXSIG=FXSI0/3. 

GO  TO  7000 

5557  X5-(A+X2)/2. 

nbr=nbr+i 
call  UN IF ( RN f NU) 


»  XSIG ) 

TO  9376 


♦  XSIG ) 

TO  9501 


t  XSIG ) 

TO  5556 


» XSIG ) 

TO  5554 


tXSTG) 

TO  5557 


,XSIG) 

TO  5553 


i XSIG) 

TO  5559 
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123 

122 


125 

124 


126 

121 


129 


123 


FOrX=XNCDF( X5,XMU, 
IF(RN.L£e  F0FX)  GO 
X6= ( X2+X4) /2. 
NBR=NRRfl 
CALL  UN  IF  <  PM , NH) 
FOFX=XNCDF{X6,XMtl, 
1 F  { R M  .  L E  .  F OF  X )  GO 
X7=( X4+X15/2. 

K!BP  =  NRP.4l 

CALL  UN  IF  (RNtA'U) 

F0FX=XNCDC(X7,XM|J, 

IFfPN.LF.FOPX )  GO 

XR=2.*X1-X7 

NRR=NRR+1 

CALL  UN  IF {RN, MU) 

FOFX=XNCDF ( X8  t  XwUt 

IFtPN.LE.FOFX)  GO 

FXMU=XR+< X8-X1 ) / 4, 

EXSIG=6.*( X3-X1  ) 

GO  TO  7000 

EXMU=( X4+X7) /2. 

EXSI G= i . 3* ( X7~X4 ) 

GO  TO  7000 

X«={ X6  +  X4) /2. 

NBR=N6R+1 

CALL  UN  IF (RM, Nil) 

FOFX=XNCOFf  XB,XM!J, 
IF(RN.LE.FQFX)  GO 
X9= ( X4+X7 ) /2  € 
NBR=NBR<-1 
CALL  UNIF(RN,NIJ1 
FOFX  =  XNCOF{X9,XM(lf 
IF(RN.LE.FOFX)  GO 
EXMU= { X4+X9 )  /  2 , 
EXSIG= 1.3*( X9-X4) 
GO  TO  7000 
EXMU=( X8+X4) /2. 
EXSIG={  X4~X8)/S. 

GO  TO  7000 
X9= ( X2  +X6  J /2 • 
N8R=NBR+1 
CALL  UNIC( PNtNU) 
FQFX=XNCDF(X<J,XMU, 
IF(RN.LE.FGFX)  GO 
EXMU=<  X6+X8)  tZ. 
EXS!G={ X3-X6J/6. 

GO  TO  700C 
0XMIJ=(  X9f  X6)  12  - 
EXSIC-=l.?*U6-X9) 

GO  TO  7000 

X7=(  X5>-X2)  /2, 

NRIUNBR+l 

CALL  UNIF(RN,MU) 

F0FX=XNC0F(X7, X^U, 

IF(RN.LF.FOFX)  GO 

XR-{  X6+X2) 

N8R=N8R+1 
CALL  UNIF(RNfN«li 
FOFX  =  XNCDF { XB  »  XMU* 
IFIRN.LEoFOFX)  GO 
X9=(X6  +  X4)/2, 
NBR=NBR  >1 
CALL  UN  IF { PN ♦NO ) 
FOFX  =  XNCOF  {  X9,  XMIJ, 
IF(RN.LE.FOPX)  GO 
EXMU=( X6+X9) /2, 
EXSIG=1 ,3*{ X9-X6) 

60  TO  7000  ^ 

EXMU={  X8*X6^X?, 
EXSIG=(X6--X8)>4. 

GO  TO  7000 
X9={  X74-X2)  (2. 


XSIG ) 

TO  5561 


XSIG  ) 
TO  121 


XSIG) 
TO  122 


XSIG) 
TO  123 


XSIG) 
TO  124 


XSTG ) 
TO  125 


XSIG) 
TO  126 


XSIG) 
TO  127 


XSIG) 
TO  128 


XSIG) 
TO  129 
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130 

127 

5561 

556  0 

9501 


9CS0 

9504 


NBR=NRP  +  1 
CAU.  UNIF(RN.NU) 
F0FX=XNCnF( X9,XMUtXSI0 ) 
IFJPN.LE.FOFX)  GO  TO  130 
E  XHtJ=  {  X2*XRJ /2. 
EXS!G=(X8-X2)/6. 

GO  TP  7000 
EXMU=( X9  +  X2 ) /2  o 
£XSIG= 1. 3*{ X2-X°) 

GO  TO  7000 
EXMU={ X7+X25/2. 

EX$IG=  1.3*U2-X7S 

GO  TO  7000 

X6=2.*A-X5 

NBR=NBR+1 

CALL  UNIF ( PNtNU) 

F0FX=XNCDF(X6,XMl)y  XSIG) 

1F(RN.LE# FOF X }  GO  TO  5560 

XR=X5 

DEL  X=X5-A 

EXMu=XRfOELX/2. 

EXS!G=1.3*DELX 

GO  TO  7000 

XB  =  X5 

DELX=A-X6 

EXMU=XQ-0ELX/4. 

EXSIG=6*f)FlX 

GO  TO  7000 

X3={ A+X2)/2« 

NBR=NBR+1 
CALL  UN I p CRM, Ml) 
F0FX=XNCDF(X3,XMU,XSIG) 
TF(RN.LE=FOPX)  GO  TO  9503 
X4= { X1+X2) /?. 

NRR=NBR+1 

CALL  UNIF(RN*NU) 

F0FX=XNCDF(X4fXMllf  XSIG  J 

JF(RM.LE.FOFX)  GO  TO  9504 

X5=2.*X1-X4 

N9R=NBR+1 

CALL  ON  IF ( RN ,NU ) 

FOFX=  X'NCDF  (  X5  »  XMU»  X3 1 G I 
IF«RN.LE.FOFX J  GO  TO  9080 
XB  =  X5 

OELX=XB-X 1 
EXMt)=XB+0FLX/4. 
EXSIG=6*0ELX 
GC  TO  7000 
XR=X  2 

DELX=X4~X2 
FXMU=XR+DFLX72. 
f:XSIG=l.3*0ELX 
GO  T.O  7000 
X5=( X3+X21/2. 

NBR=NRR+1 
CALL  UNIF(RN,MU) 

FOFX  =  X'NCDF  {  X5  ,  X*MJ,  XSIG  ) 
IF(RN. LE.FO^X J  GO  TO  9081 
X6= {  X2+X4)  /?. 

NBR=NBR+1 

CALL  l/N  Tc  (  RN »  NO } 

F0FX=XNC0F(X6tXMUt XSIGJ 

TF(RN.LE.FnPX)  GO  TO  <?0R? 

X7=2 .*X4-X6 

NBR=NRR-*-l 

CALL  UNI c { RN tM0 1 

F0FX=XNCDF(X7,X«Ut XSIG) 

IF(RN.LE.FOFX)  GO  TO  9083 

XB=X7 

0ELX=X7-X4 
EXMU=XB  tOEI.X/4# 
FXSIG-6*DFLX 
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cn  rn  7000 

9083  XB=X2 

DELX=X6-X2 
E  XMU=  XR+DEL  X/2 . 
cXSIG=l.3*DFLX 
.  '  TO  7000 

9082  FXMU-f X5+X21/2. 

EXSIG={  X2-Xr')/6. 

GO  TO  7000 
9081  X6=(X3+A}/2* 

NRR=NRR+1 
CALL  U\'IF(RN»NU) 
FQFX=XNCDF<X6, X^U, XSTOJ 
IF(RN.LE.FnFX)  GO  TO  9084 
EXMU=  (  X3*X5  }  f2. 

EXSIG= { X5-X? } /6. 

GO  TO  70G0 
9094  X7=2.*A-X6 
NaRsNpo+l 
CALL  UNIF(RN,Nt!5 
F0FX=XMCDF{X7tXMU*X8IG) 
IF(RN.LE.FOFX)  GO  TO  9085 
XR—  X  6 

DEL  X=X3-X6 
EXMU=X8+f>ElX/2. 
EXSTG=1.3*0ELX 
GO  TO  7000 
9085  XB  =  X7 

DELX=A-X7 

EXHU=xn-nFi.X/^. 

EXSI G=  6*DELX 
GO  TO  7000 
9503  X4- ( X3  +  A )  /2. 

NBR=N8R+1 
CALL  UNIF<PN,\HJ) 
F0FX=XNC0F(X4,XMU,XSTGJ 
IF(RM.LE.FOPX)  GO  TO  9507 
X5=<  X2+X3)/2. 

NBR=NBR*l 

CALL  IJMIF(RMtNt>> 

F0FX=XNCDF(X5tXMUvXSIG! 

IF(RN.LE.FOPX)  GO  TO  9509 

X6=2  >*X2-X5 

NBR=NBR+1 

CALL  IJMTP(RMfNU) 

FOFX=XNCOF(  X6f  XMi.jt  XSIG> 

!F(RN,LE.FOFX)  GO  TO  9510 

XR=:X6 

DEL  X=  X6-X2 
FXMU=XR+DELX/fc, 

EXSI G=6*DELX 
GO  TO  7000 

9510  XR-X3 
DELX=X5-X3 
EXMU=XR+DELX/2. 

EXSI G= 1 ♦ 3*9ELX 
GO  TO  7000 

v509  X6={ X3+X4) /2. 

N8R=NBR+1 
CALL  UNIF{ RNfWIJ] 

FOFX  =  XNCDF(X4»XM!l,X5IG) 
IF(RN,LE.FOPXJ  GO  TO  9511 
FXMU={  XA«-X3J  /2. 

EXSIG= ( X3-X6) /2. 

EXSIG=  2  XSIG 
EXSI G=cXS I G/ 6o 
GO  TO  7000 

9511  EXMU=( X4fX6i /?. 

EXSIC= ( X6-X4) ?? . 
EXSIG=2»*EXSIG 
EXS1G=F XSIG/6. 

GO  TO  7000 
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9507  X5=(X4+a>/2. 

NBR=NRR+1 

CALL  UNIP{RN,NIJ) 

F0FX-XNC0F( X5,XMU,XSirO 
I  F { RN,  L E  . F0C X  5  GO  TO  9508 
EXMU=< X4*X5)/2. 

EXSI G= ( X4-X5 }  /  2 . 

EX$IG=2.*FXSTG 

EXSIG=EXSIG/6. 

GO  TO  7000 

9508  EXMU=(X5*4)/2. 

EXSIG=(X4-X5) /2. 

EXSIG~2**EXSIG 

EXSIG=EXSIG/6. 

GO  TO  7000 

7000  XB=0* 

OELX=0. 

WRITF(6,100M  F  XMU  *E  XSIG 
I F( EXS IG.LT«  0. )  EXSI G=-E  XSIG 
0PT=0 . 

DIF=0. 

XI  NC  =  0 . 

0P7=EXMU-4.*EX$IG 
IF { OPT  <.LT  .  4 1  GO  TO  1003 
XINC=(0PT-AI/EXSIG 
INC=XINC/1 
DIF=XINOINC 
IFIOIF  .GE..5)  I MC  =  T NC 4- 1 
IF{ D  IF  ,LT  .  « 5 }  INCHING 
GO  TO  1001 
1003  A=OPT 
INC=0 

1001  M= 1 0  + 1 NC  + 1 
33  I S=70-N3R 
M=IQ+I NC+1 


CONDUCT  BRUCETON  TEST 


CLEAR 


10 


ARRAYS 
DO  10  1  =  1 
X(I J=0. 

ixoi  n=o 
ixxi  n=o 
Nsm=o 

NG( I  1=0 
SUMAP.  =  0  . 
SUMBR=0. 
AR  { I  )  =  0  . 
BR  ( I  )  =  0. 
CONTINUE 


200 


LOAD  X  ARRAY 

DO  20  J  =  1 1  200 
XI  J)=A  +  {  J-D+EXSIG 
2G  CONTINUE 


CONDUCT  EXPERIMENT 

30  CALL  UN  I  PIRN. Ml!) 

FOFX=X Nt'DFI  X (  M) ,X«U,XSIG) 
IF  ( RN . G7 .Fr)PX  J  GO  TO  *-0 
1 XX( M ) =1 X  X ( M ) + 1 
M=M~1 
N=N+  1 

50  IFIN.GT.IS?  GO  Tn  60 
GO  TO  30 

40  IX0(M}=IX0(M)+1 

H=M  +  1 
N=N  + 1 
GO  TO  50 


PERFORM  BRUCETQN  ANALYSIS 
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COUNT  RESPONSES  AND  NON-RESPONSES 
60  I  SllMX=0 
ISUMO=0 
DO  14  J=lt200 
I  $UMX=  T  SUMX+ 1 XX (  J  } 

ISUMO=T  SUMQ+IXO( J) 

NS { J  J  =  0 
AR{ J ) =0» 

BR( J)=0. 

NG( J I =0 

14  CONTINUE 

DETERMINE  LESS  FRFOUFNT  EVENT  AND  LOAD  NS 
IF( ISUMX.LE.tSUMQ)  GO  TO  15 
NT=!StJMO 
IFLAG=0 
DO  21  J=1 , 200 
NS(J)=IXO( J) 

21  CONTINUE 
GO  TO  16 

15  NT=  J  SUMX 
IFLAG=1 

DO  22  J=l»  200 
NSt  J  )  =  I  XX(  J ) 

22  CONTINUE 

DFfERMINE  FIRST  AND  SECOND  MOMENTS 

16  JC0UNT=1 

17  IF(NS< JCOUNT). GT.O)  GO  TO  19 
JCOUNT = JCOIJNT  + 1 

!F( JCOUNT. GF. 200)  GO  TO  104 
GO  TU  17 

18  MC0UNT=200- JCOUNT 
DO  19  J  =  1  »  MCOUMT 
NG(  J)  =  NS<  JCniJNT+J-1  ) 

AR(J)=(J-l)*NGl J) 

SUMAR=SUMA«+AR( J) 

8R ( J l=t ( J-l )**2)*NG( J) 

SUMP R =  SUM RR  +  0R ( J ) 

19  CONTINUE 
YPRIME=X< JCOUNT) 

CALCULATE  ESTIMATES  Oc  MEAN  AND  STANDARD  DEVIATION 

IF  ( I  Ft.  AG»FQ.O)  XMUE  $T= YPRI MF+EXSI G*(  C  SUMAR/NT )+(!•/£ «) ) 
I F  f  . NOT.IFLAG.FO.O)  XMUEST  =  YPR  I  Mg+EXS  I G*  C  t  SUMAR/NT } - { 1 • 
SIGFAC  =  (  (NT*3UM8R)-(SUMAR**2)  )  /(NT**2) 

IF{ S IGFAC »GT. * 3 )  GO  tq  1000 
EMU(LC0UNT)=0. 

EOEVCL COUNT) =0. 

NOGO=NQGO+i 
GO  TO  104 

10nn  OEVE  ST  =  1 ,62*EXSI G*( SIGFAC+.029) 

LOAD  EMU  AND  EDEV 

EMU(LCOUNT)=XMUPST-XMU 
EDEVtL :OUNT)=DFVFST'XSTG 
ADDMU=ADDMU<*EM(K  LCOUNT) 

ADDSIG=ADDSIG+nDEVt  LCOUNT) 

AODMUQ=AODMUO  +  FM!Jf  LCOUNT )  v*? 

ADDSDO=ADD: DO+FDEV  f  LCOUNT) «*? 

IFf»'MUCLCOUNT)  .(  T.n.  )  GO  TO  <H 
I F( f  MU( LC  nUNT ) . FQ. 0 • )  GO  TO  92 
!«UKT  =  TMUHIM 

I  F  t  FMU ( LCOUNT)  .GT.HIMH)  HI  Mtj=FMU(LCOUNT) 

Ir(  ♦NOT.FMUaCOUNT)  .GT.HIMU)  HlMIJ=Hf  «U 
GO  TO  93 
92  Nn«U=NOMtHl 
GO  TO  93 

91  I MUL  0= I MULO+ 1 

7F( EMU{ LCOUNT ) ,LT.smlO)  SMLO=EMU{ LCOUNT) 
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IF  (  . NOT , EMU(  LCOUNT )  . L  T. $ML  0 J  SMt.O=SMLO 

93  I F ( EDE V{ LCOUNT )  .  t  T . 0.  )  00  TO  94 
IF< ED£V< LCOUNT) .FO.O. )  GO  TO  «5 
IDEVHI= IDEVHI+I 

IFlEDEV(LCPUNT) .GT.QEVH1 )  DEVH I =FDFV < LCOUNT } 

IF( .NOT.PDEV(LCOUNT) .GT.dpvHI l  DFVHI =DEVH I 
GO  TO  104 

95  NODE  V=NODEV-'r  1. 

GO  TO  104 

94  !0EVL0=IDEVL0+1 

IF  ( EOHVt LCOUNT l .  LT . DEVLO)  DEVLO=EnFV(LCOUNT) 

IF{ .NOT. EDEV( LCOUNT ) «L T.DE VL O)  DEVLO=OFVLO 
104  JC0UNT=0 
SIGFAC=0. 

XMUFST =0  . 

DEVEST=0. 

SUMAR  =0  » 

$UwBR=0. 

LC0UNT=LC0UNT4-1 

HAVE  1000  EXPERIMENTS  BEEN  CONDUCTED  ? 

ifilcount.lt. icon  go  to  103 

IF  EXPERIMENTS  COMPLETED  CALCULATE  AND  WRITE  RESULTS 


EXNOGO=NOGO 

SANA VM=AODMM/( 1000.-FXN0G0) 
SAMAVD-ADQSIG/f 1 000. -F XNOGO) 
SAMSOM=ADOMUQ/{ 9T9.-FXN0G0) 
SAMSOO= ADDSOQ / ( 99^ .-EXN0G3 ) 
IFUANY.EQ.il  GO  TO  35 
IANY=I ANY+1 
GO  TO  69 
35  STOP 
END 


SUBROUTINE  UN I F  f  R N , NU l 


SUBROUTINE  RETURNS  RANDOM  NUMBER  UNIFORM  ON  (0,1). 

REAL  MOD 
MOD=  2**31 
NR=129*NU+l 
RN=NR/MPD 

IF(RN.LT.O.O)  R  N  =-  R  \* 

Nll=NR 

RETURN 

END 


FUNCTION  XNCDF{ V,XMU,SX) 


FUNCTION  SUBPROGRAM  CALCULATES  CUMULATIVE  NORMAL. 

X  IS  AN  R.V.  KITH  MEAN*  XMU, AND  STANDARD  DEV! AT ! ON, SX . 
ARG=( V-XMUJ/SX 
XNCDF=SNCDF{ ARG) 

RETURN 

END 


FUNCTION  SNCDF(X) 


FUNCTION  SUBPROGRAM  CALCULATES  STANDARD  CUMULATIVE  NORMAL. 
DATA  TEST./ 0*0/ 

I F ( TEST . NE .0 . 0 )  GO  TO  100 
SR2=  SORT ( 2.0} 

TEST= 1 . 

X  00  SNCDF=(  lcO<-ERF(  X/SR2H /2.0 
RETURN 
END 
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